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ABSTRACT 


Formulas  for  the  pseudoinverse  of  a  compact  operator 
are  applied  to  linear  system  identification  and  scattering 
function  estimation  from  a  finite  set  of  noisy  measurements. 
The  result  is  a  nonparametric  estimator  possessing  several  de- 
«  sirable  features.  The  approach  encompasses  the  Modified  Dis¬ 

crete  Fourier  Transform  and  is  applied  herein  to  the  important 
problem  of  closely  spaced  object  resolsution  in  radar/optical 
signal  processing. 
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I.  INTRODUCTION 

Extracting  information  from  noisy  measurements  is  a 
classical  problem.  Yet,  it  is  this  very  type  of  problem  which 
forms  a  long  string  of  applications  challenging  engineers  and 
scientists  of  the  past,  the  present,  and  likely  a  long  term 
future.  One  problem  area  involves  situations  wherein  an 
unknown  function  (signal)  is  observed  after  obscuration  by  a 
system  characterized  by  a  known  impulse  response  and  additive 
random  distrubances .  With  the  ever  popular  digital  computing 
machine,  the  measurements  are  recorded  digitally.  One 
therefore  faces  the  problem  of  recovering  an  unknown  continous 
function  from  a  finite  set  of  discrete  measurements. 

Motivated  by  problems  in  radar/Optics  scattering 
function  estimation  and  signal  processing,  we  study  a  solution 
to  the  above  problem  in  this  paper.  Let  s(t)  denote  the 
transmitted  signal  of  a  radar  (or  communication)  system.  This 
signal  is  scattered  by  a  target  (or  channel)  represented  by  a 
linear  system  with  impulse  response  x(t)  which  we  shall  call 
the  "scattering  function".  The  waveform  at  the  front  end  of 
the  receiver  is  the  convolution  of  x(t)  and  s(t).  Let  the 
receiver  be  a  conventional  matched  filter  with  impulse  response 
matched  to  the  transmitted  signal,  the  receiver  output  being 
therefore  a  double  convolution  denoted  as  s(t)*(s( t)*x(t) ) .  We 
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further  assume  that  the  received  data  is  recorded  in  sampled 
form.  Given  the  sample  data  set,  y(t^) ,  i-1,...,n,  we  would 
like  to  estimate  the  target  configuration  (or  channel  property) 
represented  by  the  scattering  function  x(t).  In  the  context  of 
linear  system  identification,  the  above  problem  is  the  same  as 
the  problem  of  estimating  the  impulse  response  (or  input 
signal)  of  a  linear  system  given  the  input  signal  (or  system 
impulse  response)  and  the  discrete  samples  of  the  output  of  the 
linear  system. 

Typical  approaches  to  problems  of  this  kind  fall  into 
two  distinct  categories.  The  first  one  is  to  treat  it  as  a 
deconvolution  problem.  In  this  case  the  convolution  integral 
is  approximated  with  a  convolution  summation,  and  the 
deconvolution  procedure  is  carried  out  either  via  the  discrete 
Fourier  transform  or  directly  in  the  time  domain,  as  in  [2] . 

The  second  approach  is  to  utilize  a  parametric  model  of  the 
unknown  signal  x(t)  and  attempt  to  estimate  parameter  values 
embedded  in  the  x(t)  model.  One  example  in  the  radar/optics 
area  is  the  resolution  of  closely  spaced  objects.  In  this 
case  a  maximum  likelihood  estimator  is  used  to  estimate  target 
amplitudes  and  locations  for  an  assumed  number  of  targets 
contained  in  the  observation.  In  the  latter  part  of  this 
report  we  will  use  the  closely  spaced  object  resolution  problem 
an  an  example  of  comparison  between  our  approach  and  the 


parametric  estimator  described  above. 

The  method  to  be  described  in  this  report  is  a 
nonparametric  estimator.  It  is  derived  using  the 
pseudoinversion  concept  of  operator  theory.  Our  method  is 
similar  to  the  deconvolution  method  in  that  both  utilize  a 
nonparametric  approach  and  both  are  solving  an  integral 
equation.  Our  approach  does  not,  however,  require  the 
approximation  of  convolution  integral  by  a  summation  and 
furthermore,  the  estimator  is  a  continuous  time  function  even 
though  the  data  is  given  in  sampled  form.  Our  estimator 
therefore  possesses  the  combined  features  of  deconvolution  and 
interpolation. 

We  summarize  features  of  our  approach  below.  The 
process  of  the  radar/communicaton  problem  described  above  can 
be  viewed  as  a  composite  linear  transformation  between  a  time 
function  x(t)  of  a  Hilbert  space  H  and  sampled  data  y(ti), 
i=1,...,n,  of  a  n-dimensional  Euclidean  space  En.  This 
composite  operator  consists  of  convolution  integral  operators 
with  kernels  s(t)  and  the  sine  function.  Let  Q:  H+En  denote 
this  operator,  then  one  can  write, 

y  =  Qx  (1.1) 

where  x(t)eH  and  yeR(Q)CEn,  the  sampled  data  array.  Notice 
that  we  use  R( • )  to  denote  the  range  space  of  the  enclosed 
operator.  Letting  Q+  denote  the  pseudoinverse  of  Q,  our 
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(1.2) 


/N 

estimate  x(t)  of  x(t)  is 


The  estimate  x  possesses 


given  by 


several  desirable 


properties : 


(1)  It  has  least  norm  among  all  elements  of  H  whose 
Q-image  agrees  with  y;  hence  'x  has  an  interpolating 
spline  significance. 

(2)  When  a  prior  bound  on  ||x||  is  given,  x  minimizes  the 
maximum  error  subject  to  this  bound;  hence  'x'  has  a 
minimax  significance. 


(3)  More  generally,  if  TsH+K  is  another  operation  into  a 
Hilbert  space  K,  then  T(x)  is  the  best  (minimax) 
estimate  of  T(x),  given  the  data  y  and  a  prior  bound 
on  J  x | | .  (The  case  where  T  is  the  Fourier  transform- 
is  of  particular  interest). 


In  this  report  we  will  describe  the  theory  leading  to 
representations  of  the  pseudoinverse  Q+  and  formulas  for  'x". 


This  operator  theory  is  next  applied  in  Chapter  3  to  the 


problem  described  above.  Finally,  three  examples  are 
considered  in  Chapter  4.  The  first  of  these  is  simple 
pointwise  interpolation  from  a  finite  set  of  sample  values. 

The  resulting  formula  is  the  time-domain  version  of  the 
Modified  Discrete  Fourier  Transform  of  [3].  The  second  example 
illustrates  noise-free  reconstruction  of  a  truncated  sinusoid. 


The  last  (and  most  serious)  example  concerns  the  problem  of 
closely  spaced  object  resolution,  wherein  the  original  signal  x 
is  modeled  as  an  impulse  train.  The  performance  of  our  method 
on  noisy  data  is  discussed  and  illustrated  by  several  figures. 
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II 


REVIEW  OF  RELEVANT  OPERATOR  THEORY 


Let  Q:  be  a  compact  linear  operator  on  the 

Hilbert  space  H^  with  range  R(Q)  in  the  Hilbert  space  H2 . 
With  Q+  denoting  the  pseudoinverse  of  Q,  our  estimate  x  of 


the  solution  x  to  the  equation 


y=Q(x) 


(2.1  ) 


'x'  =  Q  +  (y) 


(2.2) 


In  section  2.1  below  we  briefly  review  general  properties  of  Q 
and  Q+.  In  section  2.2  we  focus  on  the  important  special 
case  where  Q  is  of  finite  rank;  this  includes  the  situations 
where  the  data  vector  y  consists  of  sampled  values. 


2.1  Structure  of  Compact  Operators  and  Their  Pseudoin¬ 
verses 


We  begin  by  recalling  the  basic  structure  of  a 
compact  self-adjoint  (Q=Q*)  operator. 

Theorem  2 . 1  .  (Spectral  Theorem)  Let  Q  be  a  compact  self- 
adjoint  operator  acting  on  a  Hilbert  space  H.  Then  the 
non-zero  spectrum  of  Q  has  the  form  {x^,  X^ , . . .  }  where  each  X^ 
is  a  real  eigenvalue  of  finite  multiplicity,  and  lim  X^=0  if  H 
is  of  infinite  dimension.  There  is  a  corresponding  orthonormal 
sequence  {e ^ }  of  eigenvectors  such  that 
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x  e  H 


(2.3) 


Q( x )  =5  e_.  <  x,e..  >  e.. , 

Corollary  2.1.1  (Resolution  of  the  Identity)  With  Q  as 
above,  let  be  the  orthogonal  projection  of  H  onto  span  {e ^ . 
Then 


i) 

-  1 

(strong  convergence) 

ii) 

Vk  -  Vj  - 

6.  .  P. 

*3  3 

(2.4) 

iii ) 

E  =  o 

(norm  convergence) 

This  result  shows  that  Q  can  be  synthesized  from  real 
linear  combinations  of  commuting  orthogonal  projections. 


Corollary  2.1.2  (Characterization  of  the  Range)  With  Q 
as  above,  a  necessary  and  sufficient  condition  that  y  belong  to 


R(Q)  is 

i) 

y  -1-  N(Q),  the  nullspace  of  Q,  and 

ii) 

E  X-2|<y,ej>|2  <  - 

(2.5) 

If  so,  then  all  preimages  of  y  have  the  form 


<y*e  .> 


V 


(2.6) 


for  any  u  e  N(Q) . 

Theorem  2.1  and  its  corollaries  are  standard;  their 
proofs  may  be  found,  for  example,  in  [1].  Equation  (2.6)  in 
effect  specifies  the  complete  (multivalued)  inverse  of  Q.  When 
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we  use 


x 


<y,e^> 

— rr~ 

3 


(2.7) 


as  the  estimate  for  x,  we  recover  the  component  of  x  that  is 
orthogonal  to  N(Q). 

These  results  are  not  generally  of  immediate 
applicability  because  of  the  self-adjointness  requirement.  We 
can,  however,  use  them  to  establish  analogues  of  the  expansions 
(2.3)  and  (2.6)  for  arbitrary  compact  operators. 


Theorem  2 . 2  (Singular  Value  Decomposition)  Let  Q : H 1 ^H2 
be  a  compact  operator  acting  between  the  Hilbert  spaces  H  and 

I  1 

H2.  Then  there  are  orthonormal  bases  { u^ }  of  N(Q)  and  {v^ }  of 
R(Q)  such  that 

Q(x)  =  E  °  <x,u  >  v  ,  x  e  H  (2.8) 

j  j  j  j  1 

2  2  ^ 

Hence  ta-|  r  a2  •  •  •  •  J  i-s  the  non-zero  spectrum  of  QQ  (or, 

Q*Q)  • 

The  positive  numbers  a-|  ,  a2,...  are  termed  the 
singular  values  of  Q.  Theorem  2.2  generalizes  the  familiar 
matrix  singular  value  decomposition 

A  =  VDU*  (2.9) 

where  A  is  arbitrary,  (J  and  V  are  unitary,  and  D  is  a  non¬ 
negative  diagonal  matrix.  The  decomposition  (2.9)  is  easily 
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seen  to  imply  the  polar  decomposition 

A  =  UR  (2.10) 

where  U  is  unitary  and  R  is  positive  semidef inite.  Conversely 
(2.10)  combined  with  the  Spectral  Theorem  yields  a  simple  proof 
of  (2.9),  and  (2.8)  may  similarly  be  obtained  from  Theorem  2.1 
and  the  general  polar  decomposition  for  operators. 

Corollary  2.2.1  (Characterization  of  the  Range)  With  the 
notations  of  the  preceding  theorem,  a  necessary  and  sufficient 
condition  that  y  e  R(Q)  is 


i) 

y  -L  N(Q*) ,  and 

ii) 

E  <-j-2|<y.vj>|2  <  - 

(2.11) 

Inequality  ii)  of  the  corollary  is  often  called  the  Picard 
Condition.  Vectors  y  obeying  it  constitute  the  domain  D(Q+) 
of  the  pseudoinverse  Q+  of  Q: 

D(Q+)  =  R(Q)  +  N(Q*) 

{y  e  h2:  min  | |Q(x)-y |  jexists  in  H^}. 

Any  x  minimizing  ||q(x)  -  y| |  is  called  an  extremal  solution  of 
equation  (2.1).  There  is  a  unique  extremal  solution  of  least 
norm;  this  is  the  vector 

+  ^  <YrV.> 

x  -  Q  (y)  -  £  — u.  (2.12) 

j  J 

This  last  formula  defines  the  pseudoinverse  operator  Q+  on 


its  domain  D(Q+) .  We  see  that,  unless  R(Q)  is  a  closed 
subspace  of  H2  (which,  for  compact  Q,  happens  iff  Q  is  of 
finite  rank),  Q+  is  defined  only  on  a  dense  linear  subspace 
of  H2  and  is  an  unbounded  operator.  In  this  case  equation 
(2.1)  is  ill-posed  even  if  Q  is  one-to-one  and  there  will  be  an 
inherent  ill  condition  in  any  inversion  attempt. 

Although  of  limited  paractical  value  because  of  the 
difficulty  in  knowing  all  the  singular  values  0 j ,  formula 
(2.12)  is  extremely  important  theoretically  as  a  device  for 
analyzing  the  pseudoinverse  solution.  First  of  all,  it 
demonstrates  the  precise  source  of  ill  conditioning  by  showing 
that  any  perturbation  of  y  in  the  vj  direction  is  magnified 

by  the  factor  0T1  ,  which  can  be  arbitrarily  large.  Secondly, 
the  convergence  of  numerous  iterative  and  regularization 
schemes  for  equation  (2.1)  can  be  viewed  as  particular  cases 
of  a  "filtered"  pseudoinversion  wherein  the  factors  o^-1  in 
(2.12)  are  replaced  by  new  factors  <H  0 ^ ) .  Here  the  non-nega¬ 
tive  function  <P(t)  vanishes  at  0  and  behaves  asymptotically  as 
t  1  for  t  +  +*.  In  particular  when  vanishes  on  an  interval 
[0,  6]  we  have  a  "truncated"  pseudoinversion. 

2.2  Operators  of  Finite  Rank  and  Their  Pseudoinverses 

When  measurements  of  a  signal  occur  discretely  (the 


sampled  data  case)  our  observation  operator  Q  has  range  in  En 
and  hence  is  of  finite  rank.  There  is,  of  course,  in  general  a 
fundamental  dichotomy  for  compact  operators  according  as  they 
have  closed  range  (and  hence  finite  rank)  or  not.  We  recall  in 
the  following  proposition  standard  properties  of  the  pseudo¬ 
inverse  of  any  operator  with  closed  range. 

Proposition  2.1  Let  Q:  H 1  -►H 2  be  a  bounded  linear  operator 
with  range  closed  in  H2 ,  and  let  Q+  be  its  pseudoinverse.  Then 


a) 

Q+ :  H2  ♦  Hj 

is  bounded 

b) 

(Q+)+  =  Q 

c) 

Q+QQ+  =  Q+ 

d) 

QQ+Q  =  Q 

e) 

QQ+  =  orthogonal  projection  on  R(Q) 

f) 

I-Q+Q  *  orthogonal  projection  on 

N(Q) 

It  follows 

from  part  f)  that  Q+Q 

is  the  orthogonal 

projection  of  H.  on  N(Q)  and  hence  that  Q+  maps  R(Q)  (indeed, 

1 

all  of  H2 )  onto  N(Q)  .  With  reference  to  equation  (2.1)  we  may 
say,  when  R(Q)  *  H2  (as  may  be  assumed  when  Q  is  of  finite 
rank),  that  Q+  recovers  the  component  of  x  in  N(Q)"^,  although, 
strictly  speaking,  Q+  does  not  "know"  about  x.  Equivalently, 
Q+(y)  is  the  element  of  least  norm  in  the  variety  {xeHj: 

Q(x)-y}  and,  as  such,  has  an  interpolating  spline  significance. 
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That  is,  if  "smoothness"  of  elements  of  Hi  can  be  thought  of 
as  varying  inversely  with  the  size  of  the  norm  on  Hi ,  then 
Q+(y)  is  the  smoothest  interpolant  in  Hi  of  the  data  y. 

This  remark  is  particulary  pertinent  when  Hi  is  a  reproducing 
kernel  or  functional  Hilbert  space. 

Another  interpretation  of  the  estimate  =  Q+(y)  is 

possible  when  a  prior  bound  on  ||x||  i-s  given  in  addition  to 

/\ 

equation  (2.1).  Then  it  is  not  hard  to  see  that  x  minimizes 
the  maximum  error  (as  measured  by  the  norm  on  Hi)  subject  to 
this  bound.  Precisely,  x  is  the  center  of  the  "hypercircle"  in 
Hi  defined  by  the  data  and  the  prior  bound.  Thus  'x'  has  a 
minimax  significance. 

The  next  result  leads  to  the  fundamental  formula  for 
Q+  in  the  finite  rank  case.  With  it  we  can  bypass  formula 
(2.12)  and  reduce  to  a  matrix  inversion. 

Theorem  2.3  a)  Let  Q:Hi-*-H2  be  any  bounded  operator 
with  closed  range.  Then  a) 

Q+  =  Q*(QQ*)+  (2.13) 

b)  If  R(Q)  =  H2  then  QQ*  is  invertible  and  so 

Q+  -  Q*(QQ*)_1.  (2.14) 
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Part  a)  of  the  Theorem  can  be  proved  by  first  noting 

that  N(QQ*)  =  N(Q*)  and  hence  that  R(QQ*)  =  R(Q) ,  since  R(Q)  is 
closed.  Next  both  sides  of  (2.13)  are  seen  to  annuli  N(Q*)  = 
R(Q)-k  Finally,  using  the  invertibil ity  of  Q  on  N(Q)-^-  one 
shows  that  the  right  side  of  (2.13)  agrees  with  this  inverse  on 
R(Q) .  Part  b)  is  immediate. 

It  remains  to  identify  the  operators  appearing  on  the 
right  side  of  (2.14)  when  Q  has  finite  rank.  In  this  case  we 
may  assume  that  R(Q)  =  En  where  n  =  rank(Q).  If  {e^f...fe  } 
is  the  standard  unit  vector  basis  in  En,  there  exists  elements 
{qj,...,q  }  (not  necessarily  orthonormal)  in  H1  such  that 

n 

Q(x)  =  E<*,q,>e,,  x  e  H  (2.15) 

j  =  1  J  J 

Note  that  (2.15)  is  not  generally  the  singular  value 
decomposition  of  (2.8),  but  rather  is  chosen  for  convenience  in 
problem  formulation.  If,  in  fact,  the  expansion  (2.8)  is  known 
then,  as  we  have  seen,  Q+  is  directly  given  by  (2.12). 

Proposition  2.2  With  Q  given  in  the  form  (2.15)  we  have 

★  n 

Q  (y)  =  Z  <y»e .>q. , 
j  =  1  3  J 


and  the  matrix  of 


[<qj  i>  1 » 


(2.16) 


* 

QQ  = 


relative  to  the  basis  {ej,...en}. 

It  follows  that  for  y  =  (yif-»yn)T  e  En, 


Q+(y ) 


with 


(2.17) 


[<qj  r  q± >1 *  *  y.  ( 2 . 1 8 ) 

The  matrix  appearing  in  (2.16)  and  (2.18)  is  the  transposed 
Gram  matrix  of  the  vectors  {q^,...,qn}  and  is  exactly  the  Gram 
matrix  when  the  underlying  scalars  are  real.  As  such,  this 
matrix  is  positive  definite,  although  possibly  ill-conditioned, 
the  latter  depending  on  the  relative  lengths  and  alignments  of 
the  vectors  qj  in  Hi . 

When  ill-conditioning  is  present  and  there  is  noise 

in  the  data  vector  y,  recourse  may  be  had  to  a  regularization 

of  the  pseudoinverse.  In  this  case  the  coefficients  Aj  in 

0  0 

(2.17)  are  replaced  by  A^  where  e>0  and  A  satisfies 

(QQ*  +  e  I) Ae  =  y.  (2.19) 

The  corresponding  solution  in  Hi  then  has  a  smoothing  spline 
interpretation;  it  varies  in  N(Q)-^  between  Q+(y)  (for  e+0) 
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and  the  zero  vector  (for  £■*■*).  This  procedure  is  equivalent  to 
a  particular  choice  of  the  filtering  functions  4»  that  were 
introduced  at  the  end  of  section  2.1  in  the  context  of  pseuo- 
inversion  via  singular  value  decomposition.  Namely, 

<Mx)  =  -5 — . 
x  +e 


Some  examples  of  the  use  of  such  regularized  pseudoinversions 
are  given  in  section  4.3  below. 


III.  THE  PROBLEM  OF  LINEAR  SYSTEM  IDENTIFICATION/SCATTER- 

ING  FUNCTION  ESTIMATION 

We  now  apply  the  results  of  Section  2  to  the  problem 
of  linear  system  identification.  We  will  still  use  the 
scattering  function  estimation  problem  as  a  guideline  in 
formulating  the  composite  operator.  The  linear  system 
identification  problem  is  clearly  contained  in  this  problem 
area. 

3.1  The  Composite  Operator 


As  discussed  before,  the  output  of  a  radar  receiver 
can  be  written  in  terms  of  a  double  integral. 


y ( t )  - y*dx2s(xrx2)x(x2) 


(3.1) 


where  x(*)  is  the  scattering  function  to  be  estimated  and  s(») 

* 

is  the  transmitted  signal  .  The  output  y(t)  is  then  sampled  to 
obtain  the  discrete  data  y(t^),  i=1...,n.  We  can  represent  the 
sampling  process  using  a  reproducing  kernel  (RK)  K(t,x)  where 


K( t, x)  =  sinc( fl( t- x ) ) 
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The  bandwidth,  0,  is  dictated  by  the  bandwidth  of  y(t)  (not 
x(t)).  Notice  that  K(t,t)  is  a  RK  because  sine  (  •)  is  a 
function  of  positive  type  (its  Fourier  transform  is 
nonnegative);  hence 

y ( t )  =  <K(t,  •),  y( *)>  (3.3) 


whenever  y  is  3-bandlimited.  In  particular. 


yi  =  y(ti)  =y  K(t.,T)y(T)dt. 


(3.4) 


Assuming  that  s(t)  is  symmetric  about  its  mean  and  letting 
P(t)  denote  the  autocorrelation  function  of  s(t),  one  then  has 


P(t)  =  J s(t-t)s(  x)  dt 

Combining  Eqs.  (3.2),  (3.4),  and  (3.5),  one  obtains 

•x)x( x)d  x 


’i  '/p<t i" 


for  i=1 , . . . ,n. 

Applying  equation  (2.15)  to  (3.6),  one  obtains 


<  ,  P,> 


<  ,  P  > 

n  . 


(3.5) 


(3.6) 


(3.7) 


where  P^*p(t^-x)  and  the  inner  product  is  specified  by  (3.6) 
It  is  now  clear  that  Q  is  a  mapping  from  to  En. 
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Given  y^,  i=1,...,n,  to  find  x(t)  is  a  problem  of  mapping 

n  O 

from  E  to  L  .  Equation  (3.6)  also  clearly  indicates  that  this 
is  the  same  problem  as  the  linear  system  (or  input) 
identification  (or  estimation)  problem. 

3.2  The  Adjoint  Operator 

Using  proposition  2.2,  one  obtains 
*  n 

Q  (y)  =  £  y.-  p.-  (3.8) 

i  =  1  1  1 

or 

Q*  =  (3.9) 

Notice  that  the  right  side  of  equation  (3.7)  is  a  column  vector 

where  each  component  is  a  functional  represented  by  an  inner 

product.  The  adjoint  Q*  maps,  however,  from  sampled  data 
n  2 

space  E  back  to  a  space  L  of  functions  of  continuous  time. 
Equations  (3.8)  and  (3.9)  clearly  reflect  this  fact  since  pj_ 
is  a  function  of  continuous  time. 

3.3  The  Pseudoinverse  Operator 

Applying  the  above  results  to  Eq.  (2.14),  the  pseudo¬ 
inverse  operator  readily  follows: 

Q+  =  Q*(QQ*) 

=  IP1 . pnn<pi,pj>]"1  (3.10) 
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where  <P^,Pj>  is  the  (i,j)th  component  of  the  n  by  n  matrix  of 
(QQ* )  and 

<p£f  pj>  =y*p(ti-T)  p(tj -T)dT  (3.11) 

In  the  problem  of  radar  scattering  function  estimation,  P(t)  is 

the  autocorrelation  function  of  the  transmitted  signal,  and  is 

therefore  symmetric  about  t=0.  In  the  context  of  linear  system 

identification  and  input  estimation,  P(t)  is  either  the  system 

inpulse  response  or  the  input  signal  used  to  probe  the  unknown 

system,  and  is  usually  a  nonsymmetric  function  of  time.  In  any 

*  -1 

event,  the  matrix  of  (QQ  )  can  be  pre-computed  and  its 
function  is  simply  to  weight  the  data  vector  Let  a  denote 

the  weighted  data  vector,  i.e., 

a  =  (GO*)"1}:  ;  (3.12) 

then  the  final  estimate  of  x  takes  the  form 

A  + 

X ( T )  =  Q^ 

n 

55  E  P(tri)a.  (3.13) 

i  =  1  1  1 

/s 

Notice  again  that  the  estimate  x(t)  is  a  function  of  continuous 
time  t. 

The  conventional  approach  to  the  problem  of  decon- 
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volution  is  to  first  approximate  the  convolution  integral  with 
a  summation.  In  our  problem,  although  the  data  is  recorded  in 
discrete  form,  the  estimate  can  be  evaluated  at  any  instant  of 
time.  Our  method  includes  therefore  both  deconvolution  and 
interpolation.  We  will  apply  the  interpolation  significance  of 
our  approach  to  present  an  interpolation  formula  for  finite 
number  of  samples  in  the  example  section. 

3.4  Discussion  of  Numerical  Problems 

Since  the  computation  of  Q+  involves  a  matrix 
inversion,  numerical  problems  may  become  an  issue.  This  can  be 
made  clear  by  examining  Eq.  (2.12).  When  the  condition  number 

(  =  the  ratio  of  the  largest  to  the  smallest  eigenvalue)  of 
QQ*  approaches  the  dynamic  range  of  the  digital  computer  used 
for  implementation,  its  inverse  can  no  longer  be  represented 
exactly.  In  addition,  even  before  the  condition  number 
approaches  the  range  of  a  computer,  a  slight  disturbance  in  y 
(e.g.,  measurement  noise)  will  be  amplified  by  the  small 
eigenvalues  (see  Eq.  (2.12)).  Methods  for  reducing  such  effects 
are  therefore  of  practical  importance. 

Among  these  methods  are  the  truncated  and  regularized 
pseudoinverse  already  mentioned  in  Section  2,  where  it  was  also 
noted  that  both  methods  are  special  cases  of  a  filtered 
pseudoinversion  obtained  by  replacing  aT1  in  (2.12)  by  4>(  a.  ) . 


19 


Each  of  these  methods  involves  a  trade-off  between  fidelity  to 


the  data  and  noise  suppression.  As  we  insist  on  greater 
fidelity  we  amplify  the  noise  transmitted  to  the  solution,  and 
conversely.  The  final  choice  of  trade-off  may  be  made  either 
on  the  basis  of  prior  knowleldge  of  noise  statistics  or 
interactively. 

For  example,  the  truncated  singular  value  estimate 


xfi  is  obtained  from  (2.12)  by  deleting  those  terms  in  the  sum 


for  which  o.  <  5;  thus 
3 


x  6 


o^6 

3~ 


<y,v.> 

a . 

3 


(3.14) 


This  is  equivalent  to  the  ordinary  singular  value  pseudoinverse 


applied  to  the  projection  of  the  data  onto  span  {v^:  6}. 


A  key  point  here  is  that  the  method  is  independent  of  the  data 
y,  depending  only  on  the  singular  system  of  the  operator  Q. 


Hence  if  the  noise  has  large  components  along  the  remaining  o_. 


directions  in  (3.14),  there  is  unavoidable  trouble.  When  it 


can  be  assumed  that  y  includes  an  additive  white  noise 


perturbation  with  known  variance  s  ,  then  a  common  approach  is 


to  chose  the  largest  <5  that  makes  the  residual  | | Q ( x 5 ) — y | |  < 
[rank  (Q)  ]  •  s. 

A 

The  regularized  estimates  x£  were  defined  through 


Eqs.  (2.17)  and  (2.19): 


It  was  noted  that  while  xe  could  be  viewed  in  terms  of 
filtering  the  singular  value  expression  (2.12)  for  Q+(y),  for 
computational  purposes  this  expansion  is  not  required  to  be 
known. 

The  estimate  xe  can  alternatively  be  viewed  as  the 
minimizer  of  the  function 

| |Q(x)  -  y| | 2  +  e| | x | |2.  (3.16) 

In  particular,  when  the  measurement  vector  is  corrupted  with 
additive  noise  with  covariance  R,  a  weighted  minimum  norm 
solution  is  sought,  i.e., 

I  |Q(X)  "  y|  |  V1  +  G|  |x|  |2  (3.17) 

£ 

This  results  in  the  following  expression  for  X  , 

L(L-1QQ*L"T  +  ei)LTXe  =  y  (3.18) 

T 

where  LL  A  R,  as  opposed  to  (2.19).  In  this  case,  (3.15)  and 
(3.18)  give  the  regularized  estimate^. 

Regularization  is  closely  related  to  the  "ridge 
regression"  method  of  statistics,  except  that  the  rank 
assumptions  on  the  observation  operator  Q  are  somewhat 
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different  in  that  context.  Intuitively,  one  replaces  an 

initial  ill-conditioned  inversion  problem  Q(x)=y  by  a 

neighboring  well-conditioned  problem  Qe(x)  =  y,  solves  this  for 
*  *  + 

x  and  then  verfies  lim  x  =  Q  (y)  as  e  +  0.  As  such,  regulan- 
zation  is  one  type  of  representation  of  the  pseudoinverse 
Q+.  This  latter  topic  has  been  extensively  studied  for 
general  bounded  operators  [4]  . 

A  good  discussion  of  the  relative  merits  of  truncated 
and  regularized  pseudoinversion,  particularly  with  respect  to 
the  noise  level  in  the  data  is  given  in  [5], 

In  the  next  section,  three  examples  will  be 
presented.  The  last  example  deals  with  the  closely  spaced 
object  resolution  problem  where  noisy  measurements  will  be 
used.  Numerical  problems  mentioned  above  will  be  further 
examined,  and  the  results  of  applying  the  regularization  method 
will  be  studied.  The  criterion  for  choosing  the  value  of  the 
regularization  coefficient  is  in  terms  of  a  desired  false  alarm 
probability. 


IV 


EXAMPLES 


Three  examples  are  presented  in  this  section.  In  the 
first  example,  we  apply  the  theory  of  Section  3  to  derive  an 
interpolation  formula  for  finite  number  of  samples,  i.e.,  given 

A 

z(t^),  i=1,...,n,  find  z(t}  for  t  e  (-a,b)  where  a  and  b  can  be 
arbitrarily  large.  In  the  second  example,  we  use  a  truncated 
sine  function  as  the  scattering  function  and  demonstrate  the 
reconstructed  function.  The  third  example  deals  with  the 
problem  of  resolving  closely  spaced  targets  (signals). 

4.1  A  Formula  for  Interpolation  With  Finite  Number  of 
Samples 

Let  z(t)  be  a  bandlimited  function  defined  as 
z ( t )  e  BL(fl) 

=  fz(t) eL2(-",«):  Z ( w)=0  (a.e.)  for  | w | > n }  (4.1) 

where  Z(w)  is  the  Fourier  transform  of  z(t). 

From  Eq.  (3.3)  the  sampling  operation  is  represented 
by  the  reproducing  kernel  K(t,x)  given  in  (3.2),  that  is, 

K(t,x)  -  sinc(«(t-T))  ,  (4.2) 

and  the  samples  of  z(t)  are  obtained  as 


From  (3.10)  we  have  the  estimate 


z(t)  =  Q* (QQ* )_1z  (4.4) 

T 

where  >z=[z(t^)/.../z(t)] 

Q*  =  [K(tr  x)#...#K(tnf  x)J 

QQ*  =  [<K(t.,T),  K(tj,T)>] 

=  [ K ( t i  ft  j  )  1  . 

Suppose  that  the  samples  are  spaced  using  the 
"Nyquist"  rate,  i.e.. 


then 


ti+1-ti  =  ^  =  a  constant, 
<K(ti,x),  K(t  j,  t)>  =  6.jr 


(4.5) 


where  is  the  Kronecker  delta  function.  We  therefore  have 

QQ*  a  diagonal  matrix,  or 


a  n 

Z(T)  =  £ 

i*1 


z(tA  ) 


sinc( n(ti-x) ) 


(4.6) 


If  the  samples  are,  however,  taken  other  than  (faster 
or  slower)  the  Nyquist  rate,  QQ*  is  no  longer  a  diagonal 


matrix 
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.  Let  r^j  denote  the  (i,j)th  entry  of  (QQ  )  and 
a  =  (QQ* )_1z,  then 

ai  *  k?,  rikz(tk>  (4-7) 

Hence,  the  interpolated  z(t)  is 

a  n 

z(t)  =  2]  a.  sin(  fi(t  .-t)  ) 
i  =  1  1 

=  Z  E  r.kz(tk)sin(n(t.-T))  (4.8) 

i=1  k=1 

Equations  (4.6)  and  (4.8)  give  an  interpolation  formula  for 
finite  number  of  samples.  We  make  the  following  remarks: 

1)  When  sampled  at  Nyquist  rate,  the  interpolation 
formula  (eq.  (4.6))  is  the  traditional  sine  function 
weighted  with  samples.  Each  data  sample  appears  only 
once. 

2)  When  sampled  at  other  than  the  Nyquist  rate,  the 
samples  are  first  summed  with  weightings  determined 
by  (QQ*)-1.  The  sine  function  at  a  time  instant 
carries  several  data  samples. 

3)  The  above  reconstruction  formula  (with  or  without 
noise  corruption)  is  optimum  in  the  minimax  sense 
(see  section  2). 


4)  The  numerical  difficulties  of  this  approach,  if  any, 
lie  in  the  inversion  of  QQ*.  If  the  condition 
number  of  QQ*  is  too  large,  an  anomolous  response 
begins  to  occur  in  the  reconstructed  function 
(especially  with  noisy  measurements),  and  one  must 
resort  to  such  methods  as  truncated  or  regularized 
pseudoinversion  as  indicated  in  Section  3.5. 

5)  The  results  above  have  apparently  been  previously 
discussed  in  the  literature  in  the  context  of  the 
Modified  Discrete  Fourier  Transform  (MDFT) ,  [3].  Our 
results  are,  however,  obtained  with  an  entirely 
different  and  more  general  operator-theoretic 
approach.  In  particular,  our  solution  is  obtained  in 
the  time  domain,  while  that  of  [3]  is  obtained  as  a 
constrained  optimization  in  the  frequency  domain. 

4.2  A  Truncated  Sine  Function 

Let  x(t)  be  a  truncated  sine  function  over  two 
cycles,  i.e.. 


x(t)  =  sin  2*t  ;  0£t£2 

=  0  ;  otherwise.  (4.9) 

We  choose  a  Gaussian  shaped  function  to  represent  the  signal 
autocorrelation  function,  i.e.. 


P(t) 


-f2 


tv  2 


(4.10) 
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There  are  two  reasons  for  using  (4.10).  First  it  is  simple  to 
manipulate  analytically.  Second,  it  closely  approximates  the 
Hamming  weighted  and  compressed  linear  frequency  modulated 
waveform. 

We  will  only  examine  the  noise-free  case  for  this 
example.  Results  are  shown  in  Figs.  4.1  to  4.3  for  a  total  of 
9,  17,  and  25  samples,  respectively.  Notice  that  measurements 
are  taken  over  the  time  interval  where  x(t)  is  not  zero.  The 
solid  curve  gives  the  true  function.  The  diamond  shaped 
symbols  represent  the  discrete  measurements,  y^ ,  i=1,...,n. 

The  small  squares  are  the  estimates.  Notice  that  the  higher 
the  number  of  samples  is,  the  closer  the  estimates  are  to  the 
true  curve.  With  only  9  points  over  a  two  cycle  period,  the 
estimates  are  still  very  close  to  the  truth.  Since  the 
measurements  are  obtained  as  the  convolution  of  x(t)  and  p(t), 
they  do  carry  information  about  x(t)  outside  the  measurement 
interval.  This  is  evident  via  the  fact  that  the  estimates 
rapidly  go  to  zero  as  the  true  function  does. 

4.3  The  Closely  Space  Object  Resolution  Problem 

We  next  consider  a  problem  of  practical  interest  in 
radar  and/or  optical  signal  processing.  In  this  example,  x(t) 
becomes  a  series  of  impulse  function  with  each  impulse 
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representing  a  point  target  and  the  magnitude  of  the  impulse 
giving  the  signal  amplitude. 

The  Fourier  transform  of  an  impulse  function  has  a 
constant  amplitude  over  all  frequencies.  The  output  y(t)  may 
however,  be  bandlimited  as  determined  by  the  signal 
autocorrelation  function  P(t).  We  use  a  Gaussian  shaped 
function  to  represent  P(t)  as  shown  in  Eq.  (4.10)  Although  a 
Gaussian  shaped  function  is  not  bandl imitted ,  we  use  a  +  3 
standard  deviation  width  of  its  spectrum  to  approximate  its 
"bandwidth".  Using  (4.10),  one  can  quickly  obtain  the  (pseudo) 
Nyquist  sampling  interval  as  At  =  1. 

In  Fig.  4.4,  we  present  the  noise-free  measurements  and 
the  estimate  of  a  point  target.  Notice  that  the  true  function 
is  an  impulse  function  represented  by  a  vertical  line  in  Fig. 
4.4.  The  estimated  waveform  resembles  a  sine  function.  A 
circle  is  drawn  at  the  peak  of  the  recovered  pulse.  Recall 
that  an  impulse  can  be  expressed  as  the  limit  of  a  sine 
function  with  its  bandwidth  approaching  infinity.  The  estimate 
can  be  made  closer  to  an  impulse  if  more  and  denser  samples  are 
taken.  This  gradually  introduces  numerical  problems,  however, 
since  the  condition  number  of  QQ*  will  become  large. 

In  Fig.  4.5,  we  present  the  results  for  two  targets 
separated  by  1.5  standard  deviations.  For  radar  (optical) 
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TIME 


Point  target  case  (At=1) 


signals,  this  corresponds  to  a  3/4  resolution  (detector  width) 
separation.  In  this  noise  free  case,  the  estimated  function 
clearly  indicates  the  presence  of  two  targets.  Notice  that  the 
two  signal  peaks  corresponding  to  the  two  target  location 
estimates  are  indicated  by  two  circles  enclosing  these  peaks. 

In  Fig.  4.6  we  add  noise  samples  to  the 
measurements.  The  noise  samples  are  meant  to  represent  the 
white  additive  noise  at  a  radar  receiver.  At  the  matched 
filter  output  the  noise  samples  become  correlated,  with 
correlation  function  being  the  same  as  P(t).  A  signal-to-noise 
ratio  (SNR)  of  20  is  used.  This  is  an  amplitude  ratio,  i.e., 
when  the  noise  variance  is  unity,  the  signal  amplitude  is 
simply  SNR.  Results  shown  in  Fig.  4.6  indicates  that  large 
sidelobes  result  from  noisy  measurements.  This  is  due  to  the 
inherent  form  of  the  pseudoinversion  shown  in  (2.12),  where  a 
small  singular  value  ( o j )  can  cause  a  small  variation  in  y  to 
become  a  large  variation  in  x. 

In  order  to  reduce  the  noise  effect,  we  apply  the 
regularized  pseudoinversion  filter.  Notice  that  in  this  case 
we  do  not  arrive  at  an  exact  solution  for  y  =  Qx;  instead,  we 
obtain  a  solution  which  also  pays  attention  to  the  smoothness 
of  x  via  the  penalty  term  using  the  norm  of  x.  Because  our 
measurements  are  corrupted  with  correlation  noise,  the 
regularized  estimate  is  obtained  via  (3.15)  and  (3.18). 


The  results  are  given  in  Fig.  4.7.  Notice  that  the  regulariza¬ 
tion  factor  can  substantially  reduce  the  noise  effect.  In  this 
case,  the  e  is  choosen  to  be  a  tenth  of  the  inverse  signal-to- 
noise  ratio.  For  SNR=20,  e=.005.  With  the  regularization 
factor  included,  the  pseudoinversion  filter  also  includes  the 
normalization  using  the  noise  covariance  matrix. 

Also  shown  in  Fig.  4.7  is  a  algorithm  for  identifying 
scattering  centers  (targets)  imbedded  in  the  measurements. 

This  algorithm  is  a  "peak-picker"  operating  over  the  interval 
where  the  measurements  are  significantly  above  noise  (i.e.,  for 
estimates  higher  than  the  two  horizontal  lines  at  the  outside 
of  the  mainlobe).  All  peaks  identified  within  that  interval 
are  further  screened  based  on  their  amplitudes.  Only  those 
peaks  within  a  certain  range  of  the  largest  peak  are  retained 
(i.e.,  the  lower  horizontal  line  in  the  center).  This  is  done 
to  reduce  sensitivities  to  the  sidelobes  of  the  estimates.  The 
number  of  peaks  identified  gives  the  estimated  number  of 
targets.  The  peak  locations  then  become  the  target  location 
estimates.  Because  the  recovered  signal  is  a  continuous 
function  of  time,  the  target  location  estimate  can  be  nearly 
made  free  of  sampling  granularity  errors. 

We  use  Figs.  4.8-4.11  to  further  study  the  effects  of 
sampling  interval  and  noisy  measurements.  Recall  that  results 
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Fig.  4.7.  Two-target  with  noisy  measurements,  with 
regularization. 
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Fig.  4.10.  Two-target  case,  noisy  data  with  twice  the  Nyquist 
rate  and  1/2  resolution  all  separation. 
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Fig.  4.11.  The  same  conditions  as  Fig.  4.10  except  with 
regularization  in  the  filter. 
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of  Figs.  4.4-4. 7  used  the  "Nyquist"  sampling  interval  (=1). 

This  interval  is  reduced  to  half  for  results  in  Figs. 

4.8-4.11.  We  have  also  enlarged  the  horizontal  scale  so  that 
the  interval  between  samples  remains  the  same  for  the  figures. 
Noise-free  data  were  used  for  Figs.  4.8  and  4.9.  Results  of 
Fig.  4.8  should  be  compared  with  those  of  Fig.  4.5.  Notice 
that  the  resolution  is  substantially  improved  with  the 
increased  sampling  rate.  In  Fig.  4.9,  the  target  separation  is 
reduced  to  a  1/2  resolution  cell,  yet  two  targets  can  still  be 
clearly  seen.  Noisy  samples  with  a  SNR  of  20  and  .75  resolu¬ 
tion  cell  separation  are  used  in  Fig.  4.10.  It  is  evident  that 
the  pseudoinversion  filter  now  acts  as  a  noise  amplifier.  This 
is  because  when  the  samples  are  taken  closer,  the  condition 
number  of  QQ*  also  increases,  and  the  result  becomes  more 
suseptible  to  measurement  perturbations.  In  Fig.  4.11,  a 
regularization  factor  of  .1/SNR  was  applied  to  the  same  set  of 
data  used  in  Fig.  4.10.  The  improvement  is  dramatic. 

We  now  study  the  performance  of  the  closely  spaced 
objects  resolution  problem  statistically.  The  probabilities  of 
correct  indentif ication  of  target  numbers  in  terms  of  target 
separation  with  target  SNR  as  parameter  are  shown  in  Fig. 

4.12.  These  results  are  obained  with  100  Monte  Carlo  . 
repetitions.  A  regularization  parameter  of  e  =  .1/SNR  was  used 
for  all  these  results.  This  parameter  was  chosen  to  insure  a 


Probability  of  CSO  Detection 


Probability  of  Correct  Identification 


Fig.  4.12.  Probability  of  correctly  identifying  two  targets  in 
a  cluster. 
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less  than  1%  false  alarm  probability  for  a  point  target.  We 
remark  that  although  techniques  for  realtime  selection  of  e 
based  upon  the  residuals  and/or  the  norm  of  x  are  available  and 
implemented  in  many  scientific  subroutine  packages,  we  believe 
that  the  above  choice  is  application  dependent  and  a 
pre-selected  coefficient  can  substantially  reduce  computational 
burden.  Comparing  these  with  results  of  [6]  where  a  maximum 
likelihood  estimator  in  conjunction  with  the  Akaiki  criterion 
was  used,  it  is  found  that  our  method  gives  poorer 
performances.  This  is  somewhat  expected  bacause  the  method  of 

[6]  is  parametric  and  near  optimum  (except  the  use  of  Akaiki 
information).  Although  our  estimate  is  optimum,  the 
identification  procedure  is  rather  arbitrary  and  furthermore  we 
do  not  assume  a  parametric  model  for  the  targets.  Our  results 
are,  however,  substantially  better  than  the  conventional 
rule-of-thumb  resolution  criterion,  i.e.,  one  resolution  cell 
separation.  Furthermore,  our  method  is  much  easier  to 
implement  compared  with  that  of  [6] .  The  target  location 
estimation  accuracy  as  a  function  of  separation  to  a  nearby 
interfering  target  for  a  two-equal-ampltude-target  model  with 
SNR  of  20  is  shown  in  Fig.  4.13. 

Also  traced  is  the  analytical  model  given  in  [7] .  The  model  of 

[7]  is  based  upon  the  Cramer-Rao  bound  analysis  for  unbiased 
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Estimation  Error  (resolution  cell) 


Object  Separation  (resolution  cell) 

Fig.  4.13.  Location  estimation  accuracy. 
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estimates  together  with  a  bias  model  for  the  case  when  targets 
are  unresolved.  The  total  error  is  the  root-sum-square  of 
these  two  errors.  The  result  of  Fig.  4.13  is  somewhat 
suprising  because  our  error  is  slightly  smaller  than  that 
indicated  by  the  Cramer-Rao  bound  over  a  certain  region.  This 
is  because  the  Cramer-Rao  bound  is  only  applicable  to  unbiased 
estimates  while  our  estimates  are  slightly  biased  (away  from 
each  other)  thus  resulting  in  smaller  random  errors. 

In  conclusion,  we  believe  that  the  method  described 
above  is  very  attractive  for  realtime  implementation  because  of 
its  simplicity  and  good  performance.  Although  the 
identification  probability  is  poorer  than  that  of  a  near 
optimum  algorithm,  its  simplicity  in  implementation  should  be  a 
factor  for  trade-off  considerations.  Furthermore,  the 
identification  probability  is  still  within  the  bound  of 
performance  usually  considered  for  realtime  systems  [7] .  The 
estimation  error  of  our  method  is  comparable  and  even  slightly 
smaller  than  that  predicted  by  the  Cramer-Rao  bound  due  to  the 
existance  of  a  slight  bias  error. 


V 


CONCLUDING  REMARKS 


In  this  report,  we  applied  the  pseudoinversion 
operator  to  the  estimation  of  linear  systems  (or  input 
functions)  with  discrete  measurements.  Relevant  operator 
theory  was  first  reviewed.  Theorems  were  stated  with  outlines 
of,  or  references  cited  for,  their  proofs.  The  solution  to  the 
problem  of  scattering  function  estimation  was  explicitly 
given.  Three  examples  were  presented  to  illustrate  our 
results.  In  one  example,  an  interpolation  formula  for  finite 
number  of  samples  was  given.  This  result  is  the  time  domain 
counterpart  of  the  Modified  Discrete  Fourier  Transform  known  in 
the  literature.  In  another  example,  the  problem  of  closely 
spaced  targets  (signals)  resolution  was  studied. 

Since  the  pseudoinversion  involves  matrix  inversion, 
numerical  problems  may  exist  for  some  problems.  Methods  such 
as  truncated  or  filtered  singular  value  decomposition 
and  regularization  can  be  applied.  General  techniques  for 
selecting  a  truncation  or  regularization  are  available,  while 
otners  may  be  very  specific  to  particular  applications. 
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